Regresión
Jorge de la Vega
4 de octubre de 2018
Ejemplo de regresión completo
En este lab veremos qué es lo que tenemos que hacer para ajustar un modelo de regresión lineal múltiple y llegar a conclusiones finales
Datos: Accidentes de tráfico
Consideremos los datos Highway. Estos datos corresponden a un artículo de maestría no publicada de Carl Hoffstedt. Se relacionan a la tasa de accidentes de auto (en accidentes por millón de millas de vehiculos) y su asociación con varios términos potenciales. Los datos incluyen 39 secciones de Highways americanos grandes en el estado de Minnesoya en 1973. El objetivo original de estudio fue entender el impacto de variables de diseño: acpts, slim, Sig y shld que estaban bajo control del Departamento de Highways, en los accidentes.
Las variables en los datos son:
adt= conteo de tráfico diario promedio en milestrks= volumen de camiones como porcentaje del volumen totallane= número total de lineas de tráficoacpt= número de puntos de acceso por millasigs= número de intercambios señalizados por millaitg= número de tipos de intercambio de caminos por millaslim= límite de velocidad en 1973len= longitud de segmentos de Highway en millaslwid= ancho de lines, en piesshld= ancho del acotamiento en pieshtype= Una variable indicadora del tipo de camino o la fuente de fondeo del camino:- “mc” para colectores mayores,
- “fai”para highways interestatales Federales
- “pa” para highways de arteria principal
- “ma” para highways arteriales mayores
- “mc” para colectores mayores,
ratetasa de accidentes de 1973 por millón de millas de vehículos
head(Highway)## adt trks lane acpt sigs itg slim len lwid shld htype rate
## 1 69 8 8 4.6 0.00 1.20 55 4.99 12 10 fai 4.58
## 2 73 8 4 4.4 0.00 1.43 60 16.11 12 10 fai 2.86
## 3 49 10 4 4.7 0.00 1.54 60 9.75 12 10 fai 3.02
## 4 61 13 6 3.8 0.00 0.94 65 10.65 12 10 fai 2.29
## 5 28 12 4 2.2 0.00 0.65 70 20.01 12 10 fai 1.61
## 6 30 6 4 24.8 1.84 0.34 55 5.97 12 10 pa 6.87str(Highway)## 'data.frame': 39 obs. of 12 variables:
## $ adt : int 69 73 49 61 28 30 46 25 43 23 ...
## $ trks : int 8 8 10 13 12 6 8 9 12 7 ...
## $ lane : int 8 4 4 6 4 4 4 4 4 4 ...
## $ acpt : num 4.6 4.4 4.7 3.8 2.2 24.8 11 18.5 7.5 8.2 ...
## $ sigs : num 0 0 0 0 0 1.84 0.7 0.38 1.39 1.21 ...
## $ itg : num 1.2 1.43 1.54 0.94 0.65 0.34 0.47 0.38 0.95 0.12 ...
## $ slim : int 55 60 60 65 70 55 55 55 50 50 ...
## $ len : num 4.99 16.11 9.75 10.65 20.01 ...
## $ lwid : int 12 12 12 12 12 12 12 12 12 12 ...
## $ shld : int 10 10 10 10 10 10 8 10 4 5 ...
## $ htype: Factor w/ 4 levels "mc","fai","pa",..: 2 2 2 2 2 3 3 3 3 3 ...
## $ rate : num 4.58 2.86 3.02 2.29 1.61 6.87 3.85 6.12 3.29 5.88 ...Paso 1: ve los datos
Lo primero que hay que hacer es examinar la gráfica de dispersión de puntos, y seleccionar algunas transformaciones a los predictores originales con la intención de (1) lograr normalidad de los datos
hw <- Highway %>% select(-htype)
pairs(hw,col=Highway$htype) # No se grafica la variable categórica "htype"; se usa para marcar los datos
summary(Highway)## adt trks lane acpt
## Min. : 1.00 Min. : 6.000 Min. :2.000 Min. : 2.20
## 1st Qu.: 5.00 1st Qu.: 8.000 1st Qu.:2.000 1st Qu.: 6.95
## Median :13.00 Median : 9.000 Median :2.000 Median :10.30
## Mean :19.62 Mean : 9.333 Mean :3.128 Mean :12.16
## 3rd Qu.:24.00 3rd Qu.:11.000 3rd Qu.:4.000 3rd Qu.:14.60
## Max. :73.00 Max. :15.000 Max. :8.000 Max. :53.00
## sigs itg slim len
## Min. :0.0000 Min. :0.0000 Min. :40 Min. : 2.960
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:50 1st Qu.: 7.995
## Median :0.0900 Median :0.1300 Median :55 Median :11.390
## Mean :0.4005 Mean :0.2964 Mean :55 Mean :12.884
## 3rd Qu.:0.5800 3rd Qu.:0.3600 3rd Qu.:60 3rd Qu.:17.800
## Max. :2.5100 Max. :1.5400 Max. :70 Max. :40.090
## lwid shld htype rate
## Min. :10.00 Min. : 1.000 mc : 2 Min. :1.610
## 1st Qu.:12.00 1st Qu.: 4.000 fai: 5 1st Qu.:2.630
## Median :12.00 Median : 8.000 pa :19 Median :3.050
## Mean :11.95 Mean : 6.872 ma :13 Mean :3.933
## 3rd Qu.:12.00 3rd Qu.:10.000 3rd Qu.:4.595
## Max. :13.00 Max. :10.000 Max. :9.230Algunas observaciones:
La variable
sigstiene muchos 0s, que corresponden a los caminos que no tienen intercambios señalizados. Entonces se puede sustituir la variable transformando por \((sigs1 = sigs*len+1)/len\) para que tome valores positivos que pueden transformarse con potenciasPodemos ver que las variables
lwidylane, aun sin ser categóricas, toman sólo algunos valores (la primera sólo valores en 11-13 y la segunda valores 2-4-6-8). Con poca variabilidad, pueden explicar muy poco de la respuesta. Podemos ya sea considerarlas factores, o removerlas del análisisSe mecionó en clase que las variables \(v\) que tienen \(max(v)/min(v) \geq 10\) se pueden transformar a logaritmos para escalar los números. En particular
adtacptylencumplen la condición.Podemos ver en el último renglón de la matriz, que hay algunas relaciones marginales modestas con la respuesta.
Muchos de los predictores tienen algún tipo de asociación entre ellos.
Para atender las observaciones anteriores, hacemos las siguientes transformaciones y volvemos a graficar
hw <- hw %>% mutate(sigs1 = (sigs*len+1)/len,
log_adt = log(adt),
log_acpt = log(acpt),
log_len = log(len)) %>%
select(-sigs,-acpt,-len,-lwid,-lane,-adt)
hw <- hw[,c(5,1:4,6:9)] #reordena las variables en el dataframe
pairs(hw,col=Highway$htype)
Paso 2: Aplica regresión “gráfica”
Les voy a dar las ideas básicas desde un punto de vista aplicado: Supongan que la verdadera relación entre la variable de respuesta y los predictores a través de la siguiente relación: \[ E(y|X) = g(\beta'X) \] para alguna función \(g\) (completamente desconocida y no especificada). Si esto es cierto y \(y\) depende de \(X\) sólo a través de la combinación lineal de \(y\) vs \(\beta'X\), entonces la gráfica de \(\{y,\beta'X \}\) nos puede dar un indicio de quién es \(g\). Las condiciones para que este resultado (Li y Duan 1989) funcione son las siguientes:
Si los predictores están linealmente relacionados, o aproximadamente la relación entre los predictores es lineal en el scatterplot, o mejor aun, son mnormales, entonces podemos usar la regresión como la combinación lineal que necesitamos.
Dado lo anterior, las transformaciones de los predictores deben buscar linealizar ó normalizar la relación entre predictores. Esta tarea puede ser laboriosa sin herramientas interactivas, porque hay que hacer muchas gráficas.
Por ejemplo, en el scatterplot de los datos que tenemos, particularmente las variable itg y sigs1 al relacionarse con las otras variables muestra patrones no lineales. Podemos intentar transformando para ver si mejora la relación lineal:
hw <- hw %>% mutate(itg2 = sqrt(itg), #experimentando llegué a esta tarnsformación
log_sigs1 = log(sigs1),
htype = Highway$htype) %>%
select(-itg,-sigs1)
pairs(hw, col=Highway$htype)
Paso 3: Ajusta el modelo de regresión y grafica \(y\) vs \(\hat{y}\).
En este paso podemos ver cuál sería la función \(g\) a escoger, lo que sugeriría la transformación de la variable de respesta. Si la gráfica muestra una clara tendencia no lineal, entonces la respuesta debe ser transformada para corresponder a esa relación no lineal.
m1 <- lm(rate ~ ., data=hw) # incluye al factor que habíamos sacado. Estos factores usualmente no se transforman
summary(m1)##
## Call:
## lm(formula = rate ~ ., data = hw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.91012 -0.62102 -0.02195 0.69901 1.75345
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.307788 4.061093 2.538 0.0172 *
## trks -0.062485 0.092180 -0.678 0.5036
## slim -0.062432 0.061986 -1.007 0.3228
## shld -0.006614 0.108084 -0.061 0.9517
## log_adt -0.603072 0.406737 -1.483 0.1497
## log_acpt 1.064998 0.521566 2.042 0.0510 .
## log_len -0.914522 0.366275 -2.497 0.0189 *
## itg2 0.767352 0.998767 0.768 0.4490
## log_sigs1 0.714586 0.261965 2.728 0.0111 *
## htypefai 1.214930 1.589018 0.765 0.4512
## htypepa -0.856320 1.058392 -0.809 0.4255
## htypema -0.200021 0.862725 -0.232 0.8184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.028 on 27 degrees of freedom
## Multiple R-squared: 0.8095, Adjusted R-squared: 0.7319
## F-statistic: 10.43 on 11 and 27 DF, p-value: 4.178e-07plot(hw$rate,fitted(m1))
lw1 <- loess(hw$rate ~ fitted(m1),span = 0.15)
j <- order(hw$rate)
lines(hw$rate[j],lw1$fitted[j],col="red")
La transformación de la respuesta puede ser sugerida por la transformación de Box-Cox. Podemos obtener la potencia \(\lambda\) del siguiente modo:
a <-boxCox(m1)
a$x[which(a$y==max(a$y))]## [1] -0.02020202Entonces la transformación sugerida para la respuesta es \(\lambdahat \approx -0.02\) y esta transformacion incluye en su intervalo de confianza al logaritmo, así que podemos transformar la respuesta y volver a ajustar el modelo:
hw <- hw %>% mutate(log_rate = log(rate),
htype = Highway$htype)
m2 <- lm(log_rate ~ . -rate, data=hw)
summary(m2)##
## Call:
## lm(formula = log_rate ~ . - rate, data = hw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48644 -0.10418 0.00738 0.12390 0.37630
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.424825 0.969455 3.533 0.0015 **
## trks -0.021560 0.022005 -0.980 0.3359
## slim -0.019220 0.014797 -1.299 0.2050
## shld -0.007655 0.025802 -0.297 0.7690
## log_adt -0.163968 0.097095 -1.689 0.1028
## log_acpt 0.184751 0.124507 1.484 0.1494
## log_len -0.228973 0.087436 -2.619 0.0143 *
## itg2 0.115995 0.238423 0.487 0.6305
## log_sigs1 0.172316 0.062536 2.755 0.0104 *
## htypefai 0.306300 0.379327 0.807 0.4264
## htypepa -0.218802 0.252657 -0.866 0.3941
## htypema -0.148325 0.205948 -0.720 0.4776
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2455 on 27 degrees of freedom
## Multiple R-squared: 0.8002, Adjusted R-squared: 0.7189
## F-statistic: 9.833 on 11 and 27 DF, p-value: 7.58e-07plot(hw$log_rate, fitted(m2))
lw2 <- loess(hw$log_rate ~ fitted(m2),span = 0.15)
j <- order(hw$log_rate)
lines(hw$log_rate[j],lw2$fitted[j],col="red")
Este modelo parece ser suficientemente adecuado para ajustar, aunque el coeficiente de determinación es menor que el del modelo m1. Previo a la última transformación, quizá sería adecuado analizar si no hay outliers y puntos influenciales que estén alterando el ajuste. Podemos analizar los residuales como siguiente paso, en el modelo m1.
residualPlots(m1)
## Test stat Pr(>|Test stat|)
## trks 1.6715 0.10663
## slim 1.0171 0.31850
## shld -0.1145 0.90970
## log_adt -0.1793 0.85909
## log_acpt 1.7525 0.09148 .
## log_len 1.1042 0.27961
## itg2 -1.0653 0.29656
## log_sigs1 1.1115 0.27652
## htype
## Tukey test 1.4371 0.15069
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1residualPlots(m2)
## Test stat Pr(>|Test stat|)
## trks 0.4448 0.6601
## slim -0.7055 0.4868
## shld -0.6013 0.5528
## log_adt 1.1134 0.2757
## log_acpt 0.0399 0.9685
## log_len 0.3596 0.7220
## itg2 -0.1911 0.8499
## log_sigs1 0.3652 0.7179
## htype
## Tukey test -0.8103 0.4178Una variación de las gráficas de residuales, son las gráficas de modelos marginales. Esta función hace una prueba de especificación: agrega un término cuadratico a la relación entre los predictores y los residuales para verificar si la relación es curva o no.
residualPlots(m1)
## Test stat Pr(>|Test stat|)
## trks 1.6715 0.10663
## slim 1.0171 0.31850
## shld -0.1145 0.90970
## log_adt -0.1793 0.85909
## log_acpt 1.7525 0.09148 .
## log_len 1.1042 0.27961
## itg2 -1.0653 0.29656
## log_sigs1 1.1115 0.27652
## htype
## Tukey test 1.4371 0.15069
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1residualPlots(m2)
## Test stat Pr(>|Test stat|)
## trks 0.4448 0.6601
## slim -0.7055 0.4868
## shld -0.6013 0.5528
## log_adt 1.1134 0.2757
## log_acpt 0.0399 0.9685
## log_len 0.3596 0.7220
## itg2 -0.1911 0.8499
## log_sigs1 0.3652 0.7179
## htype
## Tukey test -0.8103 0.4178También podemos qq-plots, que están adaptados al caso de RLM
qqPlot(m1,id.n=3) #identifica las tres observaciones con los residuales más grandes
## [1] 5 26qqPlot(m2,id.n=3)
## [1] 34 36Por último, podemos anañizar los posibles puntos infuenciales:
influenceIndexPlot(m2)
influencePlot(m2) #combina los residuales estudentizados con distancias de Cook y valores de la matriz sombrero.
## StudRes Hat CookD
## 25 -1.7252819 0.3423982 0.12034379
## 34 -2.5622986 0.2785552 0.17514351
## 36 -1.9192166 0.1934702 0.06697469
## 38 0.7499658 0.5857165 0.06735770
## 39 -0.7499658 0.5857165 0.06735770Al parecer podemos intentar rehacer el análisis eliminando la observación 34, y posiblemente tengamos un mejor ajuste
hw34 <- hw[-34,]
m2a <- lm(log_rate ~ . - rate , data=hw34)
summary(m2a)##
## Call:
## lm(formula = log_rate ~ . - rate, data = hw34)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42397 -0.08137 -0.01081 0.13773 0.37102
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.546381 0.884011 4.012 0.000453 ***
## trks -0.006385 0.020893 -0.306 0.762333
## slim -0.025276 0.013679 -1.848 0.076048 .
## shld 0.016735 0.025349 0.660 0.514947
## log_adt -0.167700 0.088422 -1.897 0.069046 .
## log_acpt 0.153026 0.114044 1.342 0.191259
## log_len -0.198476 0.080500 -2.466 0.020598 *
## itg2 0.101466 0.217171 0.467 0.644237
## log_sigs1 0.200212 0.057973 3.454 0.001908 **
## htypefai 0.227898 0.346749 0.657 0.516802
## htypepa -0.300035 0.232231 -1.292 0.207734
## htypema -0.075962 0.189640 -0.401 0.692018
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2235 on 26 degrees of freedom
## Multiple R-squared: 0.8315, Adjusted R-squared: 0.7602
## F-statistic: 11.67 on 11 and 26 DF, p-value: 1.886e-07plot(hw34$log_rate, fitted(m2a))
boxCox(m2a)
En este último modelo podemos ver que algunos predictores no son significativos, por lo que podemos evaluar eliminarlos, considerando el ajuste de todos los modelos con estos predictores para obtener el mejor.
step(m2a)## Start: AIC=-104.29
## log_rate ~ (rate + trks + slim + shld + log_adt + log_acpt +
## log_len + itg2 + log_sigs1 + htype) - rate
##
## Df Sum of Sq RSS AIC
## - trks 1 0.00467 1.3036 -106.154
## - itg2 1 0.01091 1.3098 -105.973
## - shld 1 0.02177 1.3207 -105.659
## <none> 1.2989 -104.291
## - log_acpt 1 0.08995 1.3888 -103.746
## - slim 1 0.17056 1.4695 -101.602
## - log_adt 1 0.17970 1.4786 -101.367
## - log_len 1 0.30369 1.6026 -98.307
## - htype 3 0.52532 1.8242 -97.384
## - log_sigs1 1 0.59584 1.8947 -91.943
##
## Step: AIC=-106.15
## log_rate ~ slim + shld + log_adt + log_acpt + log_len + itg2 +
## log_sigs1 + htype
##
## Df Sum of Sq RSS AIC
## - itg2 1 0.01444 1.3180 -107.736
## - shld 1 0.03091 1.3345 -107.264
## <none> 1.3036 -106.154
## - log_acpt 1 0.08757 1.3911 -105.684
## - slim 1 0.18064 1.4842 -103.223
## - log_adt 1 0.18514 1.4887 -103.108
## - log_len 1 0.32341 1.6270 -99.733
## - htype 3 0.54355 1.8471 -98.910
## - log_sigs1 1 0.68196 1.9855 -92.165
##
## Step: AIC=-107.74
## log_rate ~ slim + shld + log_adt + log_acpt + log_len + log_sigs1 +
## htype
##
## Df Sum of Sq RSS AIC
## - shld 1 0.02881 1.3468 -108.914
## <none> 1.3180 -107.736
## - log_acpt 1 0.07856 1.3966 -107.536
## - log_adt 1 0.17997 1.4980 -104.872
## - slim 1 0.19074 1.5087 -104.600
## - log_len 1 0.35547 1.6735 -100.662
## - htype 3 0.84669 2.1647 -94.882
## - log_sigs1 1 0.72374 2.0417 -93.104
##
## Step: AIC=-108.91
## log_rate ~ slim + log_adt + log_acpt + log_len + log_sigs1 +
## htype
##
## Df Sum of Sq RSS AIC
## <none> 1.3468 -108.914
## - log_acpt 1 0.12455 1.4714 -107.553
## - log_adt 1 0.16303 1.5098 -106.572
## - slim 1 0.18728 1.5341 -105.967
## - log_len 1 0.42026 1.7671 -100.594
## - htype 3 0.83505 2.1819 -96.582
## - log_sigs1 1 0.70358 2.0504 -94.943##
## Call:
## lm(formula = log_rate ~ slim + log_adt + log_acpt + log_len +
## log_sigs1 + htype, data = hw34)
##
## Coefficients:
## (Intercept) slim log_adt log_acpt log_len
## 3.27719 -0.01962 -0.13859 0.16805 -0.22216
## log_sigs1 htypefai htypepa htypema
## 0.20563 0.29093 -0.28367 -0.09170El modelo sugerido es el que tiene el AIC más pequeño.